Chiral d + is superconducting state in the two dimensional t-t' Hubbard model 
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Applying the recently developed variational approach to Kohn-Luttinger superconductivity to the 
t-t Hubbard model in two dimensions, we have found, for sizeable next-nearest neighbor hopping, 
an electron density controlled quantum phase transition between a d-wave superconducting state 
close to half filling and an s-wave superconductor at lower electron density. The transition occurs 
via an intermediate time reversal breaking d + is superconducting phase, which is characterized by 
nonvanishing chirality and density-current correlation. Our results suggest the possibility of a bulk 
time reversal symmetry breaking state in overdoped cuprates. 



As demonstrated by Kohn and Luttinger (KL), the 
generic ground state of a clean degenerate electron sys- 
tem should be superconducting^ Direct application of 
the KL approach to the t-t' Hubbard model in two di- 
mensions has shown that, at sufficiently weak coupling 
when no particle-hole instabilities develop, the phase di- 
agram of the model in the t'/t versus electron density, 
p, plane, exhibits several superconducting phases of dif- 
ferent symmetry^ Similar conclusions have been reached 
previously in mean field studies of the t-t 1 -J model^ and 
more recently making use of the renormalization group 
methodic and fluctuation exchange approximation. 7 
Thus, by tuning one control parameter (t'/t or p), a 
quantum phase transition between two different super- 
conducting states can be realized. The purpose of this 
paper is to study such transitions within a simple but 
realistic microscopic model. 

As a particular example motivated by the physics 
of high temperature superconductors, we have chosen 
to study the two-dimensional t-t' Hubbard model with 
the bare single particle dispersion = — 2t(cos k x a + 
cos k y a) + At' cos k x a cos k y a where a is the in-plane lat- 
tice constant. This is the minimal model including 
electron correlations and reproducing the shape of the 
Fermi surface observed by angle resolved photoemission 
spectroscopy^ in the cuprates. In a previous study^ we 
have shown that close to half filling d-wave is the domi- 
nant pairing symmetry in a broad range of t'/t and that, 
under reducing the electron concentration, different pair- 
ing symmetries start to dominate. Here we have cho- 
sen to study the model for t'/t = 0.45, in which case 
a density-controlled transition between the d-wave and 
s-wave pairing states was predicted to occur in the vicin- 
ity of p w 0.7i£ This choice of t'/t is motivated by the 
observation^ that at smaller ratios of t'/t, the transition 
either happens at small electron fillings, or is located 
close to the Van Hove density in which case our method 
does not apply. Moreover, according to a recent study, 9 
Tl and Hg based materials do indeed support sizeable 
next nearest neighbor hopping. 

We address this problem within the recently devel- 
oped variational approach to superconductivity. 10 The 
main idea of Ref. [l(J is as follows: Performing a canon- 
ical transformation which eliminates the scattering pro- 



cesses to first order in U/t (where U is the local Hubbard 
interaction), we have constructed an effective Hamilto- 
nian which leads to an interaction in the Cooper chan- 
nel of the form Vk P = U + C/ 2 xi(k + p, e p £k), where 
Xi(q, <jS) is the real part of the particle-hole susceptibil- 
ity x(q, w) = L- 1 2 k (/k - /K+q)/(£K+ q -ejc-u-iO) 
and L is the number of lattice sites. Superconductivity is 
realized if a nontrivial order parameter Ak which solves 
the gap equation 



k kp p 2E„ 



(1) 



can be found. In Eq. (|TJ we introduced the quasiparticle 
enegy E p = (£ 2 + lApl 2 ) 1 / 2 , where £ p = e p - p and p is 
the chemical potential. 

All calculations reported in this paper were performed 
for moderate interaction strength U/t — 7. Strictly 
speaking, this is outside the region of validity of our 
weak-coupling approach. The need for such values of U 
is due to the finite system sizes, typically up to 512 x 512, 
which we can study within a reasonable time. At smaller 
values of U, the finite size effects become appreciable. 
However, we do not believe the large values of U are a 
serious drawback of our calculations. In fact, in the d- 
wave state the value of U is just a multiplicative factor 
entering Vk P (since the d-wave state is not affected at 
all by the constant term U in Vk P ). On the other hand, 
although both terms in V^p do affect the s-wave pair- 
ing state, we have checked a posteriori that our solutions 
in the s-wave channel exhibit very small on-site pairing 
amplitudes, thus providing a nontrivial self-consistency 
check of the calculation. 

At p — 0.685 we have solved the gap equation at T = 
by an iteration method. The calculation was performed 
many times, starting from a random complex vector Ak. 
We have found three types of solutions: pure d-wave, 
pure s-wave, and d + is-type solutions which can be writ- 
ten as Ak = dk + *Sk with real dk and Sk- As a next step, 
we have required that the self-consistent solution is one 
of the above three types (which speeded up the calcula- 
tion considerably) and in this way we calculated the data 
presented in Fig. ^i. Note that pure s-wave and d-wave 
solutions are stable for all studied densities, whereas the 
d + is-type solution can be stabilized only in the interval 
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FIG. 2: Density of states per lattice site (in units of t) for 
p — 0.69 in the pure d-wave, pure s-wave, and in the d + is 
state (dashed, thin solid, and thick solid lines, respectively). 
Calculated following Ref.[l(|for the same parameters as Fig.0 
The finite weight at low energies is due to the finite width of 
the S functions, j/t = 2 x 10~ 4 . 
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FIG. 3: Magnitudes of the real (left panel) and imaginary 
(right panel) parts of the pairing function, \Dij\ and \Sij\, in 
the chiral d + is state. The middle of the square corresponds 
to i = j. Calculated for p = 0.69 and the same parameters as 
in Fig. [JJ,. 



FIG. 1: (a) Maximal value of |Ak| at the Fermi surface. 
Triangles: pure d-wave state. Squares: pure s-wave state. 
Circles: the same for |dk| and |sk| in the d + is state, (b) 
Condensation energy per lattice site (in units of t) in the pure 
d-wave (triangles), pure s-wave (squares), and in the d + is 
state (circles), (c) Phase diagram of the t-t' Hubbard model, 
(d) Spin chirality £123. The upper inset defines the lattice 
points 0,1,2,3, and 4. The lower inset shows the pattern of 
current chirality on the plaquette 0123. All data were calcu- 
lated as functions of the electron density p for t' /t = 0.45 and 
U/t — 7 on special lattices— with L = 512 x 512 (a,d) and 
L = 256 x 256 (c). (b) was calculated by extrapolating the 
results in (a) to 4096 x 4096. 



0.655 < p < 0.715. 

We have calculated the BCS condensation energy 
-Bcond = -Sk(^k - |6c|) 2 /2i?k as a function of p for 
all three types of solutions. The results are plotted in 
Fig. We find that the d + is state is the true ground 
state in the whole region of its stability. Thus, at weak 



coupling, the quantum phase transition between the s- 
wave state at low electron density and the high-density 
e?-wave state occurs via a mixed symmetry state d + is, 
in agreement with a model calculation for an isotropic 
Fermi surface and separable interactions^ Within our 
mean field like approach, both s/d + is and d/d + is 
transitions are of second order. Field theoretic studies 
suggest that for the d/d + is transition this is valid also 
beyond mean ReldiZ^ 

It is worth pointing out that the complex admixture 
increases the condensation energy substantially, up to 30 
per cent close to the crossing point of pure s-wave and d- 
wave states. Since the condensation energy is determined 
by the typical value of |Ak| close to the Fermi energy 
and, according to Fig. QJl, the maximal value of |Ak| 
is approximately the same in the mixed and pure states 
close to the crossing point, the density of states in the d+ 
is state must be substantially suppressed in the subgap 
region. This is indeed found in an explicit calculation, 
see Fig. |3 Note that a small true gap is opened in the 
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d + is state, whereas both the d-wave and the s-wave 
states are gapless^ 

In order to get further insight into the nature of the 
d + is state, let us consider the singlet pairing function 
Fij = <j(c.i- a Cj a ). In the d + is state we can write 
Fij = Dij 



iSij where = L 1 X)k(^ k /Ek) cos k ' ^ 



is the vector 
are 
in the 



Sij = L 1 X)k( s k/-E'k) cosk • Ry , and R 4 
connecting lattice sites i and j. Note that and Sij 
real and even in In Fig. we plot \Dij | and \S. 
chiral state. We find that \D^\ is peaked for j—i = (1,0), 
(2,0), and symmetry related points, whereas \Sij\ peaks 
for j — i = (3,0), (3,1), (2,2), and symmetry related 
points. This means that the d-wave and s-wave states 
do not compete for the same Cooper pairs, allowing for 
a coexistence of both states in the d + is state. 

The temperature dependence of the d- wave and s-wave 
parts of the gap in the d+is state is shown for p = 0.7 in 
Fig. 0^.. Two critical temperatures and T c d+is can be 
defined, the mixed symmetry state being realized at T < 
T c d+is , and the pure d-wave state occuring at Tj? +is < 
T < T c d . Collecting the data for different p — const cuts 
together we obtain the phase diagram Fig.^Ji. Note that 
the region of stability of the d + is state shrinks with 
increasing temperature, due to the lower entropy in the 
mixed state than in the pure states. 

The complex order parameter in the d + is state im- 
plies that time reversal symmetry is broken. Moreover, 
the point group symmetry is reduced from C^ v down to 
C 2v , i.e. the 90° rotations and reflections across the di- 
agonals, m', do not belong any more to the symmetry 
operations of the d + is state. A physical order parame- 
ter which distinguishes the two possible time-reversal re- 
lated states d+is and d—is can be constructed as an odd 
function of the local spin operator S 2 ; = 2~ 1 a a pc\ a Cip {a 
is the vector of Pauli matrices) or of the charge current 



15 



jij = -itij J2a( c icr c ja ~ c )a c i<?) defined on the bond ij 
However, because of spin rotation invariance of the d + is 
state, we have (S<) = 0, while the unbroken m reflection 
symmetry, i.e. x — > —x, forces (Jij) = 0. Therefore more 
complicated functions of Si and jij have to be considered. 

Let us start by taking the simplest order parameter 
made up solely of the spin operators. This must nec- 
essarily be a third power of Sj. The spin chirality 16 
Eijk — (Si ■ (Sj x Sfc)), where i, j, and k are lattice sites, 
is a natural candidate for such an order parameter, since 
it is spin rotation invariant. Very recently, the d + id xy 
RVB spin liquid has been characterized^ by Eijk aud in 
what follows we show that this order parameter is useful 
in our case as well, if the lattice points i, j, and k are 
chosen in such a way that a nonzero value of E^k is not 
prohibited by the remaining Ci v symmetry. 

In a singlet superconducting state, the chiral order pa- 
rameter can be evaluated as 



E. 



ijk 



where G, 



-Im [GijFj k F^ + GjkF ki F*j + G kl F tJ F* k ] ,(2) 



r) = L 1 J2kh cosk 
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pendent of a and an even function of i and j . For conve- 



FIG. 4: (a) Maximal amplitudes of the real and imaginary 
parts of Ak at the Fermi surface in the d+is state, (b) specific 
heat cv per lattice site, and (c) normalized penetration depth 
and spin susceptibility as functions of temperature. Calcu- 
lated for p — 0.7 (and the same parameters as in Fig. on 
a special lattice— with L = 256 x 256. (b) and (c) were 
calculated by extrapolating the data from (a) to a lattice 
1024 x 1024. 



nience let us introduce the following abbreviated notation 
for the neighborhood of the lattice point R (see inset to 
Fig. Elf): = R- x + y, l = R-x, 2 = R, 3 = R + y, 
and 4 = R + x, where x and y are unit lattice vectors 
in the x and y direction, repectively. Consider first the 
triangle 123. Since G\ 2 = G23, S12 = S23, D 12 = -D 2 3, 
and D13 = 0, making use of Eq. J2J one finds readily that 
E 123 = (3/2)D 12 (S 12 G 13 - S 13 G 12 ). In Fig. EJi we plot 
E123 as a function of doping. 

Let us point out that the chiral order parameter pro- 
vides also a direct measure of broken spatial symmetries. 
In fact, the nonvanishing value of E\ 23 is a direct con- 
sequence of broken m! symmetry, since in a m' symmet- 
ric system Ei 23 = E 32 i = — Ei 23 . Moreover, a clock- 
wise rotation of the lattice around the point 2 leads to 
-Ei 23 ~~ * E324, and a straightforward calculation shows 
E324 = —Eyi 3 . On the other hand, under the m reflec- 
tion we have -E123 — > E423 = E\ 23l consistent with the 
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FIG. 5: Density-current correlation function Bijk- Lattice 
point i is in the center. The arrow points from j to ft for 
Bijk > and from k to j otherwise. For clarity, only half of 
the diagonal currents is shown. Calculated from the pairing 
function in Fig. [3] 



invariance of the d + is state under the m reflection. 

Let us turn to order parameters based on current op- 
erators. The simplest order parameter of this type isi& 

Bijk = ^(ndjk +jjkni) = t jk lm(F*F ik ), 

and the equality applies to a general singlet supercon- 
ductor. For a d + is superconductor, this specializes 
to Bijk = tjk (DijSik - SijDik). An example of B^k is 
shown in Fig. EI The nonvanishing density-current corre- 
lation B^k suggests that potential scattering on imper- 
fections such as impurities or surfaces will, in general, 
cause finite currents and magnetic fields in the vicinity 
of these defects. For scattering on nonmagnetic impu- 
rities, such currents have been explicitly calculated only 
recently (for a, d + id xy superconductor)^ On the other 
hand, nonvanishing surface currents have been predicted 
already by Anderson and Morel^ In our case the surface 
currents should be most pronounced for [110]-like sur- 
faces, while for [100]-like surfaces the unbroken m sym- 
metry forbids the formation of such currents. This is 
fully consistent with the predictions of Ginzburg-Landau 
theory^ 

Next let us consider, in analogy to spin chirality, the 
third power of local current operators, Pijk = {jijjjkjki — 
jikjkjjji}- Note that, being an expectation value of a Her- 
mitian operator jijjjkjki + (jijjjkjki), the order parame- 
ter P^k is a real quantity. A nonzero value of Pijk means 
that the product of three current operators along a closed 
loop is different for current operators ordered clockwise 
and anticlockwise around the loop. Therefore we call 
Pijk the current chirality. A lengthy but straightfor- 
ward calculation shows that, in a singlet superconductor, 



Pijk — ~{i^/i)tijtjktkiEijk- Thus, up to a trivial mul- 
tiplicative constant, the spin and current chiralities are 
equal. Making use of the symmetry properties of #123 we 
find that the four possible current chiralities which can be 
defined on an elementary plaquette 0123 of the square lat- 
tice satisfy P103 = -P012 = -^032 = -Pias- The resulting 
current chirality pattern is shown in the inset of Fig.^i. 
Note the similarity of Fig.QJi with the distribution of cur- 
rents in a chiral state of a three-band model introduced 
by Varma as a description of the pseudogap phase of the 
cupratesi^i This similarity is consistent with the recent 
symmetry classification of superconducting states coex- 
isting with a chiral translationally invariant pseudogap 
phase of nonsuperconducting origini22 The crucial differ- 
ence of our calculation with respect to Ref. is that 
in our case no net currents are flowing, (jij) = 0. This 
makes it possible that a translationally symmetric chiral 
state is realized also in a one band model. 

Our weak coupling approach, which should be reliable 
sufficiently away from half filling, therefore predicts a 
bulk chiral d + is state in overdoped cuprates with size- 
able next-nearest neighbor hopping. Next we address 
the question about experimental observability of such a 
state. Upon cooling, two phase transitions should be vis- 
ible. Figure 0Jd shows that this is indeed the case but 
the anomaly of the specific heat at T c d+is is small, since 
the quasiparticles are already frozen out to a large extent 
at T^ +ls . This is expected to be a generic result which 
holds except for the case when T c d ~ T^ +ls . On the other 
hand, within BCS theory the in-plane penetration depth 
is given by 



A = 



£oda 2 h 2 c 2 L 



2 r 



where d is the interlayer distance. In Fig. 0fc we show 
explicitly that T^ +ls should be observable in the temper- 
ature dependence of A, and also of the spin susceptibility 
X- Note also the activated behavior at low temperatures 
which is consistent with the full gap. 

Direct tests of time-reversal symmetry breaking are 
more involved, since in a perfect sample (jij) = and 
therefore no spontaneous magnetic fields develop. More- 
over, since the photoemission experiments test only the 
diagonal part of the Nambu-Gorkov electron Green's 
function, the broken parity in our d + is state can not 
be directly measured by circular dichroismSi of photoe- 
mission spectra. 

Recently a spontaneous splitting of zero bias conduc- 
tance peaks in a tunneling experiment on overdoped 
cuprates has been interpreted in terms of a bulk time re- 
versal symmetry breaking superconducting state. 23 How- 
ever, this experiment measures the properties of a [110] 
oriented surface of a d-wave superconductor and in this 
geometry, time reversal symmetry may be broken in the 
vicinity of the surface even in absence of such symmetry 
breaking in the bulk^S Therefore phase sensitive exper- 
imental tests of the putative bulk chiral d + is state at 
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large doping should avoid geometries with in-plane edges. 
We propose to measure the current-phase relation I(<p) of 
c-axis Josephson junctions between a conventional super- 
conductor and overdoped cuprates. Assuming that the 
Josephson coupling is due to a finite s component of the 
cuprates, from the relative sign of the first and second 
harmonics of 1(f) one can determine the relative phase 
of the d and s components of the order parameter^ 

In conclusion, applying the recently developed varia- 
tional approach to KL superconductors^ we have shown 
that for sufficiently large next-nearest-neighbor hopping 
i'/t, the Hubbard model predicts a change of pairing sym- 
metry from d-wave to s-wave with overdoping, in agree- 
ment with a previous calculation^ The quantum phase 
transition between the two superconducting states occurs 
via an intermediate chiral d + is phase. We have shown 
that the density-current correlation and the spin and cur- 



rent chiralities serve as local measures of the broken time 
reversal and parity symmetry in the d + is state, and we 
calculated their spatial distribution. Our results suggest 
the possibility of a bulk time reversal symmetry breaking 
state in overdoped cuprates with sizeable next-nearest 
neighbor hopping. We have calculated thermodynamic 
quantities which might be useful for experimental tests 
of the putative d + is pairing state. 
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